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Aa Ss TRACT 

A computer simulation study of the effect of a 
thermalized target lattice on the ranges of irradiating 
ions has been made. The calculations are in good agree- 
ment with experimental results. 

The work is a continuation in the development of 
a computer model formulated at the USNPGS which 
takes into consideration the displacement of atoms in 
the target lattice as well as inelastic energy losses by 
the primary ion. The simulation was done for a Xenon 
ion striking the (100) face of a tungsten target. This 
thesis attempts to establish the relative importance of 


thermal and inelastic effects in the n=body model. 
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l. Introdaetion 

A history of the study of radiation damage should 
be subdivided in the experimental, analytical and simu- 
lation approach to this problem. Most of the earliest 
theoretical work was done in the analytical manner, 
but recently the high speed computer has become a 
useful tool, and the simulation treatment of radiation 
damage has emerged as an important adjunct to the 
other theoretical approaches. 
Theoretical 

An excellent and very readable review of the 
analytical theory is contained in Dienes & Vineyard 
"Radiation Effects in Solids"! which should be consulted 
by all serious students of radiation damage. In it, they 
note the work from the first efforts of Bohr“ in 1948 
on the interatomic potential, up to just about the time 
the first work was being done with a computer. 

Seitz and Koehler? produced a model which treated 
a locally excited lattice region as a temperature spike. 
More work on a particular range of temperature spikes, 


the highly excited type, was done by Brinkman4»>, 


This "temperature spike" model provides some of the 
prominent theoretical explanations of radiation effects. 

Another of the several valuable analytic models of 
displacement production was originally proposed by 
Kinchin and Pease®. They assumed that all collisions 
were binary, the lattice atoms were initially at rest, 
collisions between pairs of moving atoms were unlikely, 
and the thermal energy of the lattice is negligible. 

At this time early workers were encountering 
difficulties in explaining some of their results. The 
back diffusion of moving lattice atoms to the surface 
should be improbable in view of the mean free path of 
the incident and lattice ions. Silsbee’ proposed and 
Leibfried® developed an energy chain concept which 
allowed direct momentum transmission along close=packed 
directions of the lattice. With this concept, many of 
the puzzling results could be explained immediately. 

Gibson, Goland, Milgram and Vineyard? were able 
to verify this energy chain concept with the use of the 
digital computer. They used a copper atom in a model 


of the copper lattice and assumed an n=body, i.e. one 





in which the primary ion is interacting with several 
target atoms at one time. With Born-Mayer potential 
functions they were able to obtain data that approxi- 
mated low energy experimental results quite well. 

Erginsoy, Vineyard and Englert!0 continued work 
with this n=body model in the body centered cubic 
lattice, and showed that the primary created an inter- 
stitial at some distance from the point of impact and 
&@ vacancey at the original site. They also showed the 
threshold energy for this phenomenon was highly depen- 
dent on the direction of impact. 

The binary collision model has been investigated 
extensively by Robinson and Oenll, le, They first assumed 
a randomized model of a solid which was inadequate to 
explain the strongly penetrating component observed by 
the CHALK RIVER experimental group, and later a 
lattice structure model for the target. As mentioned 
before, a range dependance on the initial direction was 
observed. 


Harrison!3 »/4 and his students have used the n-body 


treatment to make further studies of the effect of a 
target lattice structure on the depth of penetration, 
scattering angle and recoil angles. Gay and Harrison!3 
compared the results of this n=-body treatment to the 
binary collision model. Harrison, Leeds and Gay|4 
found that the total range computed from the n-body 
model substantially exceeded that obtained from the 
binary collision model. 
Experimental 

Among the earlier workers who attempted heavy 
ion range measurements were Schmitt and Sharp), 
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Powers and Whaling’ and van Lint, Schmitt and 
Suffredini!? who measured the range of initial ions with 
energies in the kev ranges. 

Recently J.A.Davies and co-workers have obtained 
a great deal of information!$:19s£03;2l with an ingenious 
technique. They bombard the target with radioactive 
ion, then anodize a surface layer, strip it chemically and 
measure its radioactivity. By repetition of the anodizing 


and stripping, the fraction of bullets ions remaining after 


a specified range could be obtained. This process how- 


ever is limited to materials which could be anodized. In 
their work, they demonstrated that the nuclear reaction 
yields were extremely sensitive to crystal orientation in 
each of the fcc, bcc and diamond type lattice. 

Lutz and Sizmann** obtained similar results in an 
experiment with Kr on Cu. ‘The ion energy could be 
varied from 10 to 15 Kev and ranges were determined 
in the (lll), (110) and (100) directions. Again, the 
channeling effect was reported along certain preferred 


directions. 








a Study Objectives 


This report is a continuation of the USNPGS 
program. In their latest paper!+ Harrison, Leeds and 
Gay are careful to inform the reader of two inherent 
errors contained in their assumptions} 


The bullet moves in a perfect lattice, which is 
undisturbed by thermal displacement of the atoms. 
Preliminary calculations indicate that the thermal 
effects are too large to neglect, but their actual 
contribution is still unknown. °°*°° 


All "inelastic" energy losses have been neglected. 
Moving atoms are known to lose energy by a 
"friction" method at low velocities as well as by 
the more familiar: resonance effects which appear 
at higher energies. 


This thesis, in attempting to establish the relative 
importance of these two errors, makes a comparison of 
1,20 


the results to those obtained by Davies et a in an 


experiment conducted at Chalk River, Canada. 





3% Simulation Model 

The basic model is that originally developed by the 
USNPGS group over the past three years. A brief 
discription of this model follows but for exact details, 
the reader is referred to reference 13, 14 and the paper 
by weéde~-. 

The simulation model consists of a single primary 
and a bee target lattice. The primary can be fired 
into the target at any angle and can be aimed toward 
any desired point; however, in view of the objectives 
of this thesis, all runs were made with the primary 
striking the target normally. 

The target lattice initially is a 39 atom structure 
with each atom subject to a small random displacement 
which simulates the thermal displacement. As the ion 
proceeds through the target, the target is built on in 
front of the primary ion. This rebuilding of the target 
is accomplished by two "layers" of atoms at a time. 

At the same time the last two layers are stripped off 
and discarded. ‘The target lattice can be rebuilt on all 


six faces. 





The heart of the simulation program is the mechanism 
by which the primary ion interacts with the remaining 
atoms. ‘This model uses a double iteration procedure to 
determine the forces on the primary as it advances. The 
interatomic potential function and force function are of 
the exponential type, 

F = exp(A + Bx) 
where A and B are empirically determined constants and 
x igs the atomic separation. 

The double iteration procedure is best described by 

Harrison et al./3 


The unbalanced force....(on the primary ion) 
is an average force calculated by a double 
iteration procedure as follows: (l)assume 

an atom at position 1 with velocity 1; (2) 
calculate the total force on the atom as a 
result of all the other atoms in the lattice 
(this means normally only about 8-10 nearest 
atoms); (3) call the calculated force, force 
l, and use the equation of motion to move the 
atom to a temporary position, position @3 

(4) now repeat the force calculations for 
position 2, call this force 2; (5) go back to 
position 1, and use the average of force 1 and 
force & to move the atom to a new position, 
position 3. Procedure | through 5 constitutes 
one "time step".13 


This double iteration procedure in an n=-body model 





is more accurate than other models previously tried on 
the digital computer. 

The model used to describe the thermal vibration 
of the target atoms (the primary was not considered 
to carry any thermal energy but rather initially placed 
in a fixed position) consisted of the application of a 
random displacement to the atom from its "perfect 
fixed lattice" position. The maximum amplitude of this 
displacement was consistent with the assumed tempera- 
ture, and a triangular approximation to the Gaussian 
distribution was used to determine the probability of 
displacement !“, The random displacement was computed 
separately for each axis to give a proper spherical 
displacement. 


The approximation to the Gaussian is 


“tp , 
Pcp) = (27T) exp ora ) 


~ ig = sey, -(6 tos" 
ie @ le < lel < oo 


This distribution is multiplied by an amplitude, 


dependent on the assumed target temperature, which was 


an input parameter to the program, and a factor of 
//3 to make the displacement spherical since this 
correction is being made along each orthagonal axis. 


Thus 


b= (AN) @ - %) 
= (44)(1- &) 


where D is the displacement to be applied and A(T) is 
the thermal amplitude. 

The computer supplies a random number between 
0 and 1, this is equated to r/ 6 and the computed 
correction then applied to the basic lattice position. 

The subprogram which supplies the above random 
number is a standard program available in the "library" 
of the computer. It works with an original input 1} 
digit random odd number, supplied by the computer 
programer for the first run, and thereafter by the 
program itself. To assure the thermal correction may 
be applied in both the plus and minus direction, the 
computer checks the eleventh digit of the input para- 


meter, if this digit is odd, the correction is subtracted, 
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and if even, the thermal correction is added to the 
basic position. 

In the true physical crystal, each atom's displace— 
ment influences the adjacent atom's thermal displacement. 
Thus it is improbable that two adjacent atoms would be 
displaced to their maximum amount in opposite direction. 
Therefore as a displacement is computed for a particular 
lattice site, only seven tenths of it is applied to the 
correction and three tenths of the adjacent atom's 
displacement is applied to the correction. 

The amplitude of displacement is of course depen- 
dent on the assumed temperature of the target. This 
amplitude was computed in the manner of Wert and 
Thomson“4, The following derivation follows this 
‘reference with slight changes. 

Assume a central atom is surrounded by six atoms, 
one on each end of the three orthagonal axis. Displace 
the atom such that only the distance to the atoms on 
one axis vary (to the first approximation). Then the 


change in energy is 


AE = % fE(%- 4) + Ely,4q)- LEC)? 


I} 





where u = x = Xp and z is the number of nearest 
neighbor atoms. The factor of two (2) appears because 
the movement of the atom causes both the energy of the 
of the displaced atom and of the neighbor atom to change 
by an equal amount. ‘The first term is the energy of 
bond with the right neighbor and the second with the 
left neighbor. 

Since the function E(x) can be expanded about its 
minimum, the total energy change for small displacements 


is 


AE=(*%f ave le (x-%,)*= (xo) 4* 
wr CK) (PE) 


An atom bound to a definite site by the potential-— 
energy law of the above equation is a simple harmonic 
oscillator. Hence the force "ff" on each atom as a 
function of its displacement is 

~ = ~d(akeyy, = eae 4 
Thus for unit displacement,x is the force on the 


atom. 


l2 


The differential equation describing the atom's 
motion is 


m dou fyr? = XU 


where m is the mass of the atom. A solution is 

u = A, cos /X/m Z 
and therefore the frequency is 

_ = ts) J Yr 

can be estimated by assuming Hookes Law is applied 

to a unit cube of side a4),. 

Y 44, = (2 
Y = Youngs Modulus, —F =A, the force necessary to 
stretch the cube a unit distance. Since Young's Modulus 
is about loll newton/m*, o“% is approximately 25 newton/m. 
Classically, the energy of the atom is 


AE = 35 ae 


Thus at room temperature 
°o 
AE = Ss KxZ ao 


_2/ 
~ .2 */0 


Since 


cb 
| 
Y 


(2) u® 


ys J2AE 
“A 
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/ 24/48 ~jo7*' 
25 


A. jo ve 


lA ~v 06 Lattice oni te 


A slightly more refined calculation gives u =.05 


lattice units. 


Similarly at freezing temperature and melting 


temperatures, the amplitude is .025 and .08 respectively. 


lh 


desis Procedure 

This simulation study was conducted on the CDC 
1604 computer, using both FORTRAN 60 and symbolic 
computer language. 

The computer sets up a bcc lattice of atoms of 
tungsten, see figure l, and starts a Xenon ion into it. 
A brief study of the symmetry of the target lattice 
shows that the indicated impact triangle covers all 
possible points in the (100) face of tungsten. Within 
this impact triangle, twelve points were chosen, 
figure 2, as being representative of the triangle. Con- 
sideration of the physics of the problem as well as the 
economics of computer time were both instrumental in 
determining just how many points should be chosen in 
this impact triangle. 

The primary ion was initially fixed one lattice unit 
away from the target (perpendicular to the plane of 
the paper in figure 2) and given an energy varying 
from | to 100 Kev. 

Barly runs were made on position 2l, figure 2, 


since this the center of the channel between the two 
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nearest atoms on the face of the target. Simulated 
thermal effects showd be more meaningfw on these 
channel and edge of channel positions than on others 
where the primary strikes a target atom essentially 
head on and rebounds. Results showed this to be true. 

After making several total range runs, with an 
initial energy of 1 Kev and with both a thermalized and 
non=-thermalized lattice, it was determined that cutting 
off the primary ion at 50 "time steps" (section 3) 
would give a profile from which an accurate DE/DX 
cowd be determined. DE/DX is computed by dividing 
the difference of the primary ion's initial energy and its 
total energy after 50 time steps by the depth of pene- 
tration along the direction of firing. 

Runs were made to determine DE/DX vice the 
total range since otherwise the computer time would be 
prohibitive, and a comparison of DE/DX for thermalized 
vs non=-thermalized lattice is just as fruitful and accurate 
as a comparison of total ranges in assessing the effects 
of thermal energy. 

One hundred runs with initial energies of 1 Kev 
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were made at each of the twelve positions and the mean 
DE/DX and deviation from this mean were determined. 
F'rom this data, a sample of twenty five runs was 
determined to give a satisfactory mean DE/DX. There- 
after all samples were of twenty five runs; this was 
convenient from the viewpoint of actual time required 

on the computer as well as sufficiently accurate 


physically. 
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5- Results 

Figure 3 through l2 shows comparisons of the rate 
of energy dissipation at various initial energies for each 
of the several impact points. Position 25 and 29 were 
not reported as the primary ion, striking so close to a 
lattice atom, produced no statistically consistent data 
other than the fact it did rebound back out every time. 

It is apparent from these graphs that those 
impact positions on the edge of the open channels of 
the lattice were the most affected by the thermalizing 
of the target. Even these, once the initial energy was 
above about 10 Kev, were not affected noticeably by 
the thermalizing energy. 

For these impact points of most interest, data 
was also obtained for lattice temperature of about 
absolute zero and for the melting point temperature. 
This data was only obtained for initial ion energies of 
1 to 10 Kev, the most significant range for this effect. 

The deviation of each position's mean DE/DX and 
each position's mean penetration were recorded. Again, 


only if the point of impact was close to an atom of 
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the target did this deviation rise appreciably above 0.5% 
of the principal value. 

Finally the total range of the 5 and 20 Kev ions 
was obtained by a graphical integration. The reciprocal 
of the rate of energy dissipation was plotted against 
the initial ion energy, and knowing the total range at 
the low end of this curve, the total range at any other 
point could be found from the area under the curve. 

These ranges for the several points of impact 
were plotted on the impact triangle and contour lines 
for various penetrations were drawn, see figure 13 and 
lye Considering the total area of the triangle as unity, 
the area under any one of the contour lines represented 
the percent of particles not yet stopped at this range. 
This is the required data for the "Davies plot," figure 
1 Sims 

On this plot, experimental work taken from the 


1.°9 is also plotted. However, the 


paper by Davies et a 
curves are artificially matched at the one per cent point 


Since, quantitatively, the total range of the ion from 


any as yet devised computer program does not exactly 


Ig 





match experimental results. In addition, errors in the 
graphical integration, which would tend to all be towards 
the same side, would be corrected by this procedure. 

One final interesting point was noticed although 
not explained, figure 16. In the course of taking data, 
runs were made with various thermal amplitudes varying 
from the value for temperatures well below the 
theoretical absolute sero to slightly above the melting 
point temperature. The points, when plotted versus 
the rate of energy dissipation, produced a curve con- 
sistent with expectations except between the value for 
theoretical zero temperature and actual zero amplitude. 
Here there was a sudden dip for values of thermal 
amplitude between 107+ and 16>, The value then rose 
back up at a thermal amplitude of zero, i.e. non- 


thermalized target lattice. 
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6x Conclusions and Acknowledgements 

Thermalizing the target lattice does appear to 
make a definite difference in radiation damage studies 
in the lower energy range. However this effect is 
slight compared to the effect of "inelastic" energy 
loses whence the ion loses its energy by "friction" to 
the electron of the target. In any future study of 
radiation damage, the effect of thermal energy can be 
accounted for by proper choice of interatomic potential 
function and constants, which is still the central and 
most difficult problem to solve in the studies of this 
nature. 

A model such as the one used for this study 
appears to give correct qualitative results and is 
adequate for comparison to experimental results. 

In future studies of the thermal energy effect, 
at least one lesson can be learned from this study. 
Since computer running time is so expensive and will 
become more so, any extra subprograms which are 
unnecessary should be bypassed. ‘This thermal energy 


need not be applied to the target for points of impact 


al 


near the atoms of the target lattice, nor for targets 
where the initial ion energy is 10 Kev or greater. 

I would like to extend thanks to Professor D.E. 
Harrison Jr. for the help, guidance and advice offered 


during the preparation of this paper. 
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